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Abstract 

A revised iterative method based on Green function defined by quadratures 
along a single trajectory is proposed to solve the low-lying quantum wave function 
for Schroedinger equation. Specially a new expression of the perturbed energy is 
obtained, which is much simpler than the traditional one. The method is applied 
to solve the unharmonic oscillator potential. The revised iteration procedure gives 
exactly the same result as those based on the single trajectory quadrature method. 
A comparison of the revised iteration method to the old one is made using the 
example of Stark effect. The obtained results are consistent to each other after 
making power expansion. 
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1. Introduction 



Recently a single trajectory quadrature method is proposed to solve the low-lying quantum 
states of N-dimensional Schroedinger equation [1,2]. The ground state wave function in 
N-dimension can be expressed by quadratures along the single trajectory. Furthermore 
a Green function is defined along the single trajectory [2]. This makes it possible to 
develop an iterative method for obtaining the ground state wave function, starting from 
a properly chosen trial function. The method is applied to solve one dimensional double- 
well potential [3]. The convergence of the iterative solution very much depends on the 
choice of the trial function. 

However, the solution based on this iterative formula is not completely the same as 
the one obtained from the single trajectory quadrature method. In the single trajectory 
quadrature method the wave function for the ground state is proposed to be in the form 
of an exponential, as it should be, while in the original iterative solution the obtained 
correction to the wave function is in the form of a power expansion. This sometimes leads 
to unreasonable results, as will be shown in the following. 

In this work a revised version of the iterative solution is introduced, which is based 
on a different integral equation using the same set of Green functions. A new and much 
simpler expression of the energy correction is obtained. The solution is in an exponential 
form and the result is exactly the same as the one obtained based on the single trajectory 
method. 

In Section 2, a brief introduction is given about the Green function method based 
on the single trajectory quadrature. Special discussion is given for the revision of the 
iterative formula. A comparison of the original and the revised version of the iterative 
solution is shown for a simple example. The revised iterative formula is then applied 
to the unharmonic oscillator potential and a comparison of the two iteration solutions 
is made using the example of the stark effect in Section 3. Finally some discussions are 
given at the end. 
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2. Green Function and the Revised Iterative Solution 



(1). Green function and the revised iterative formula 

We discuss a particle with unit mass, moving in an N-dimensional unperturbed po- 
tential Vo(q). The ground state wave function $(q) satisfies the following Schroedinger 
equation: 

tf$(q) = £$(q), (2.1) 

where 

H = T + V (q) = ~V 2 + V (q). (2.2) 

Assume the solution of Eq.(2.1) could be expressed as 

$(q) = e' s ^\ (2.3) 

where S(q) is the trajectory based on which the Green function will be defined. From 
Eqs.(2.1)-(2.3) it is easy to derive the following equation for S(q): 

^(VS) 2 -^V 2 S-V + E = 0. (2.4) 

Introduce perturbed potential U (q) and assume 

V(q) = V (q) + U(q) (2.5) 
H = T + V(q) = -^V 2 + V(q). (2.6) 

Define another wave function ^(q) satisfying the Schroedinger equation 

Wtt(q) = £tt(q) (2.7) 

£ = E + A. (2.8) 

Let 

*(q) = e -S(q)-r(q)_ ( 2 .9) 

The equation for r and A could be derived easily [2]: 

VS- Vr + i[(Vr) 2 - V 2 r] = (U - A). (2.10) 
Consider the coordinate transformation 

Qi, <?2, <?3, ■ ■ ■ , Qn -> S, «i, a 2 , • • • , "v-i ( 2 -H) 
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with a = (cci, ct2, • • • , ctjv-i) denoting the set of iV — 1 orthogonal angular coordinates 
satisfying the condition 

VS'-Va i = 0, (2.12) 

for i — 1, 2, • • • , TV — 1. Similarly to the discussions in Ref. [2], introducing the ^-function 
in S-space: 

and define 



c = e[(vsy}-. (2.i4) 



Using 



VS-VC=1, (2.15) 
it is easy to derive the following equation 

r = G[(U - A) - i(Vr) 2 ] = (1 + CT)- l C[(U - A) - i(Vr) 2 ], (2.16) 

where the Green function G — (1 + CT) _1 C. From this equation and the expression of 

J$ 2 (qMq)e- T Wdq 

J$ 2 (q)e-(q)rfq ' 1 ' 

the iteration process could be performed in the following way: Starting from A = and 
ro = 0, 

_ J$ 2 (q)t/e-— 
" ~ J$2(q) e -^-irfq' 

r n = (I + CT^CW-A^-^Vt^) 2 ]. (2.19) 

When the N-dimensional variable q is transformed into (S, a), T = — |V 2 could be 
decomposed into two parts: 

T = T S + T a , (2.20) 

where T$ and T a consist only the differentiation to S and a, respectively. The detailed 
expression of T s and T a could be found in Appendix A [2]. Now another function could 
be defined as [2] 

D = -26e 2S %e~ 2S h s h a (2.21) 
ti a 
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and it is related to the Green function G and C in the following way [2]: 

G=(l + DT a )- 1 D = (l + CT)- 1 C. (2.22) 
Thus, the integral equation for r and A could also be expressed as 

r = G[(f/-A)-i(Vr) 2 ] = (1 + DT a )^D[{U ~ A) - \{Vrf] 

= ^(l + T^r^fJ-A)-^) 2 ]. (2.23) 

The explicit expression of r is 

r = -2 I" e 2S '^fdS' I*' e- 2S "hs„h a dS"(l +T a D)- 1 [(U - A) - i(Vr) 2 ]. (2.24) 
jo ri a Jo I 

Therefore, we have 

-^e- 25 ^^ = j\-™'hs>KdS\l +T a D)- 1 [(U - A) - i(Vr) 2 ]. (2.25) 
The left hand side of Eq.(2.25) approaches when S — > oo, so is the right hand side, i.e., 

roo 1 

/ e- 2S h s h a dS(l+T a D)- 1 [(U-A)--(Vr) 2 ]=0, (2.26) 
Jo 2 

which is correct for all a. Integrating over da = Tlfs^daii and because of 

hsh a T a = -- g (2-27) 



and 



we derive 



J dah s h a T a r = 0, (2.28) 



h s h a dadSe- 2S [(U - A) - ^(Vr) 2 ] = 0. (2.29) 



Denoting dq = hsh a dadS we reach a new expression of the perturbative energy 

A 



Jdq e- 2S [U-\{VTf 



J dq e~ 2S 

1 -(X7t) 2 ] 

(2.30) 



Jdq &[U-\{Vt) 2 ] 



Jdq $ 2 




(2.31) 



(2.32) 



Or, we can also use Eqs.(2.31) and (2.19), instead of Eqs.(2.18) and (2.19), as our iteration 
series, since the calculation of Eq.(2.31) is much simpler than Eq.(2.18). 

Comparing the above iteration with the old one derived from the equations for / = e~ T 
and A [2]: 



there are several advantages for the new one: 

1) It directly gives an exponential form for the perturbed wave function e~ T . This 
result is consistent with those obtained in the single trajectory quadrature method using 
the series expansion of {Si} and {Ei} (See Section 1 of Ref.[2]). 

2) It makes the iteration more transparent, since there are no terms appearing asUf 
in the integration. 

3) The calculation of the perturbation energy is much simpler. 

(2). Example 

This formula has been tried for several examples. In the following the result for a 
simple example will be given, with the one-dimensional potential 



A. 



fn 



■n 



1 + G(-U + A n )f 7 
/dq* 2 /n-i ' 



n— 1 i 



(2.33) 




34) 



Now we have 



V (x) = \g 2 x\ 



(2 
(2 



35) 
36) 



U(x) = Xx 
and according to Eqs.(2.1) and (2.3) 



S(x) = \gx\ 



(2 
(2 



37) 
38) 



From Eqs.(2.14) and (2.15) we have, in this example 



I (qx 



(gx) 2 
X 1 

= [ — dx . (2.39) 
J qx 
o y 

Define 

Starting from r = and A = 0, based on Eqs.(2.18) and (2.19), we have 

00 _ 2 

/ e gx Xxdx 

A l = =2|_ = o (2.41) 

/ e~ gx2 dx 



and 



The next step gives 



Ao = 



CC/ = ^x. (2.42) 



/ e-^Xxe~ x dx 



— oo 
l2 



A 2 



While further iteration gives 



(2.43) 



r2 = (i + CT)- 1 C[(C/-A 2 )-^(Vr 1 ) 2 ]=r 1 . (2.44) 

The iteration is completed here and the obtained corrections for the ground state energy 
and wave function are 

A - ~. (2-45) 

r = -x. (2.46) 
9 

This is the exact solution of the ground state for the potential of Eq.(2.34). 
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Performing the iteration based on Eqs.(2.31) and (2.32), starting again from r = 
and A = 0, we have 



D = -2 j e 9y2 dy J e'^dz. (2.47) 

-oo 

Ai = and 



Tl = DU = -x. (2.48) 
9 



The expression of A 2 is 



A 2 = 



1 e-^\\x-\^)dx 

-oo y 
oo 

j e~ gx2 dx 

— oo 

12 



\2 

(2.49) 



which is the same as the result in Eq.(2.43) and 

= DU = r 1 . (2.50) 

The iteration is completed and the exact solution is obtained after only two steps of 
iteration. There is no restriction of the parameter A to obtain this result. 

If applying the old iteration based on Eq.(2.33) for f n and A n , the obtained corrections 
for the ground state energy would be expressed as a ratio of two polynomials, while the 
wave function would be a polynomial. Although one seems to reach the correct energy 
correction of A 2 = — ttt in the second order of the iteration, the equation similar to 
Eq.(2.50) never appears between / 2 and /i, which means that the iteration would never 
complete. The f n becomes more and more complicated with larger n, so is A n . Only 
when A 2 << g 3 , these expressions are meaningful. Under this condition, after making the 
series expansion, the expression could approach the series expansion of the exact solution. 
Therefore, the simple expression of the exact solution could not be reached in the old 
iteration procedure. 
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3. Application of the Revised Iterative Formula 



In this section the revised iterative formula is applied to solve the anharmonic oscillator 
potential. The two iteration procedures are compared in details using the Stark effect as 
an example. 

3.1 Anharmonic Oscillator potential 

As an example we discuss the one-dimensional harmonic oscillator with perturbative 
potentials, applying the revised iterative formula based on the Green functions defined in 
Eqs.(2.14), (2.21) and (2.22). And the the result is compared to those obtained in ref.[2] 
using the perturbation expansion based on the same set of Green functions. 

Example 1. Consider a one-dimensional harmonic oscillator with 





and a perturbative potential 



U — ex 



(3.2) 



where p is a positive integer. The unperturbed ground-state wave function is 




(3.3) 



which satisfies Schroedinger equation 




(3.4) 



Since S — \ gx 2 , we have 




(3.5) 



The eigenstate of 




H — g 2 x 2 + e x 



is 



e 3 T which satisfies Schroedinger equation 




(3.6) 



For one- dimensional case, Eq.(2.22) becomes 




(3.7) 
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The revised iterative formula now consist of Eqs.(2.19) or (2.32) for r n and Eq.(2.31) for 
A n as following: 



/ dx e~ 9X \ ex 2p -\ (t^ 2 



Taking 



p 

m=l 



and 



(3.8) 



j dx e 9x2 

— oo 

rn = (l + CT)- 1 C[ e ^-A„-i«_ 1 ) 2 ] 

= Dlex^-K- 1 -^) 2 ]. (3.9) 



t o (x) = 0, A = 0, (3.10) 
based on the formula given in Appendix B, it is easy to obtain 

oo 

/ dx e~ gx e x 2p 

Ai = — = er lp , (3.11) 

/ dx e~ 9X ' 2 

— oo 

T\ = (1 + CT)- X C [ e x 2p - eT 1 p ] 
= D [ex 2p -eT lp ] 



= eJ2r mp x 2m (3.12) 



t[ = t J2 2mT mp x 2m -\ (3.13) 

m=l 

For the second order, we have 

oo 

/ dx e~ 9X [ e x 2p - \ [t[) 2 ] 



J dx e 9X 

— oo 

Y v v 

= eTi P - - e 2 H 4mnr mp r np ri >m+n _i. (3.14) 

m=l n=l 

The iteration could continue with r 2 , A 3 , etc. Finally this will lead to the energy correction 
A n and the wave function in the order of n: 

m n (x) = e~ s ~ Tn . (3.15) 
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Specially, for p = 2 



3 o 21 

2 



A 2 = e ^ - e 2 



{2gY 8gf 

The result is the same as those from the single trajectory quadrature method introduced 
earlier in ref.[2]. When expanding e~ T to a power series, it is consistent to the perturba- 
tion expansion results obtained in ref. [2] based on the same set of Green function. 

Example 2 Using the same one-dimensional harmonic oscillator Hamiltonian (3.1) as 
the unperturbed H, we consider now an odd power perturbation potential 

eU = ex 2p+1 . (3.16) 

With analysis parallel to that in the first example, using the formula introduced in Ap- 
pendix B, starting from 



it is easy to obtain 



and 



t o (x) = 0, A = 0, (3.17) 

Ai = , (3.18) 

T\ = (1 + CT)- l C [ e x 2p+1 ] 

= D[ex 2p+1 ] 

= eJ2l mp x 2m+1 (3.19) 

m=0 



t[ = e £ (2m + l) 7mp x 2m (3.20) 

m=0 



Similarly to Eq.(3.14) we have, for p ^ 0, 

Jdx e"^ 2 [ex 2p+l -\ (r[) 2 } 



Ao = 



/ dx e~ 9x2 

-9 e ' E E( 2m + X )( 2n + l)7mp7n P r 1;m+n . (3.21) 

^ m=0 n=0 



For the special case of p — 1, we have 

Ai = 



1 1 

g 2 3g 

A, = r- 



n = e (-^ x + ^- x 3 ) 
11 



1 



89 

T2 = n-e 2 (/-^--i-x 4 ). (3.22) 
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The result of p = is already shown as the simple example at the end of Section 2. These 
results are exactly the same as those obtained based on the single trajectory quadrature 
method. After expanding e~ Tn to a power series the result is consistent to the one based 
on the perturbation expansion using the same set of Green function in ref.[2]. 

It is worthy to mention that the formula based on the Green function D has another 
advantage for one-dimensional case. Since the expression of r n is a double integration, 
T' n could always be expressed as a single integration, i.e., no differentiation is needed. 
Comparing to the single trajectory quadrature method introduced in ref[2], this revised 
iteration method could reach the same result, however, it has much less restriction to the 
potential. Because the existence of many times of derivatives in the earlier method it 
requires the potential being continuous and very smooth for high order derivatives, while 
the revised iteration formula based on D is more flexible and could even be applied to 
discontinuous potentials. 



3.2 Stark Effect 



As a simple multi-dimensional example we now discuss the Stark effect, making de- 
tailed comparison of the two iteration procedures. Let H be the unperturbed Hamiltonian 
for a Coulomb potential: 

1 a 2 

H = — V 2 - — , (3.23) 

2 v 

with V 2 denoting the three-dimensional Laplacian, r the radius. Consider the Stark effect 
as perturbation the corresponding Hamiltonian is 

U = H + ercostf, (3.24) 

where i? is the polar angle; i.e., 

r 2 = x 2 + y 2 + z 2 and z = rcos-#. (3.25) 
The ground-state of H, $(r) = e~ s = e~ g2r satisfies 

H$(r) = E$(r) = --g 4 $(r) (3.26) 

and the ground-state wave function of Ti satisfies 

Wtf(r) = (E + A)tf(r). (3.27) 

Introducing 

tf(r) = e - s(r) - r(r) , (3.28) 

the perturbative energy A and the function r(r) introduced for the wave function in 
Eq.(3.28) could be solved iteratively based on the revised iterative formula Eqs.(2.19) or 
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(2.32) for r n and Eq.(2.31) for A n . The corresponding functions C and D used in the 
Green function G are expressed as (see Appendix C) 

r 

C = \j dr (3.29) 
g o 

r oo 

D = -2 j dr^e 2g2r J dr'r a e^ r '. (3.30) 

r 

Using the definitions of corresponding quantities in the spherical coordinate given in 
Appendix C, we can perform the following two iterations. Therevised one is 



Jdr e- 2 9 2r [ ercostf - \ (Vv^ 2 ] 
n ~ Jdr e-Wr ' 

r n = (l + CTy'C [ erC os^-A„-i (Vr n _i) 2 ] 



(3.31) 



= (l + DT^D [ erC os^-A„-i (Vr n _i) 2 ], (3.32) 
starting from A = and r (r) = and the old one is 

a Jdr e~ 2g2r ercos^/ n _i(r) 

n ~ Jdre-^r fn _ l{r) > 

/„ = l + a + C^-^^-erco^/^xCr) 

= l + ^ + DT^-^tAn-ercos^/n-^r), (3.34) 

starting from A = and /o(r) = 1- From the revised iteration formula (3.31) and (3.32) 
it is easy to obtain 

f dr e~ 2g2r \ ercos-tfl „ ,„ „ . 

1 4 

Tl = e( V r2 + (Vp r)cos ^' (3 ' 36) 

To calculate A 2 and r 2 we need 

(Vr '» 2 = e2[( (W + <2?F r + (W r> » + < <W r + (W rVs2 * 1 (3 ' 37) 

and the obtained results are 

Jdr e-^fercos^-KVri) 2 ] 2 36 
A2 " /dr e-^r " _e (W ( } 

r 2 = ^r 2 + ^±_r)costf 
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Further iteration will give 



A - £2 (W +e 1P +0(£6) ' (3 ' 40) 

The result is exactly the same as those based on the single trajectory quadrature method. 
The old iteration based on Eqs.(3.33) and (3.34) gives 

f dr e~ 2g2r ercosfi , , , 

Ai = — =-g = 0, 3.41 

/ dr e~ 2 9 r v ' 

A(r) = 1 " ^ ^ + Ocostf. (3.42) 



Further iteration gives 



jdr e- 2g2r ercosi?/i(r) _ 2 36 
jdr e-nW ~ ~ C (V) 4 ' 



A J m c ^ ^ j 1 y y _ 2 ^ 



/ 2 (r) = 1 - e(^ r 2 + -^Ip r)cos^ 

+e2[ 3WF r3(1 + 3cos2 ^ + Wf r2(1 + cos¥)] 

+ ^ 2 [(^, r 2 ) 2 + (^ r) 2 ]cos 2 ^ + 0(e 3 ). (3.44) 

Now we could see clearly, the perturbed wave function e~ T in the revised iteration is an 
exponential while f n in the old iteration is a polynomial, which could be meaningful only 
when e/g 2 << 1. It can be seen from Eqs.(3.39) and (3.44) that e~ T ' 2 and f 2 are consistent 
to each other up to the order of e 2 , after making a power expansion for e~ r2 . The energy 
correction for the 4th order from the old iteration formula is expressed as a ratio of two 
polynomials. However, after making the power series expansion for the factor e/g 2 , it 
gives also the same result as Eq.(3.40) given from the revised formula, up to the order 
of e 4 . Therefore, to compare the results of the two iteration procedures the power series 
expansion for the factor e/g 2 for both cases are needed. Up to any fixed order of e/g 2 the 
two results are consistent to each other, if the iteration reaches the corresponding order. 
Of course, to keep the expansion meaningful the condition e/g 2 << 1 should be satisfied, 
although the revised iteration itself does not need this condition, if one does not make 
the expansion. 
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3. Discussions 



In this paper a revised new iterative method based on Green function defined by quadra- 
tures along a single trajectory is proposed to solve the low-lying quantum wave function 
for Schroedinger equation. The method is applied to solve the unharmonic oscillator po- 
tential. A detailed comparison of the revised iteration method to the old one is made 
using the example of Stark effect. 

From above analysis the advantages of this revised iteration method can be seen clearly. 
Because it directly gives an exponential form for the perturbed wave function e~ T , it some- 
times gives better convergence than the old formula, like the case of unharmonic oscillator 
shown at the end of Section 2. The exponential form of the perturbed wave function also 
ensures the result of this revised iteration procedure consistent with those obtained us- 
ing the series expansion of {Si} and {E{\ based on the single trajectory quadrature(See 
Section 1 of Ref.[2]). When we look at the formula Eq.(2.30) for the calculation of the 
perturbation energy, it is completely new and much simpler than the one in the old it- 
eration procedure (see Eq.(2.17)). The expression of the perturbation energy in the old 
iteration method is similar to the one in the usual perturbation method. In the new 
iteration formula Eq.(2.30) the denominator does not change in each order of iteration. 
This makes the calculation much simpler. The improvement for the energy and the wave 
function in each order of iteration is very clear. This revised iteration procedure will be 
applied to more other potentials in our future work. 
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Appendix A 

For the convenience of applying the coordinate system {S, a} the definition introduced 
in ref. [2] is given in the following. Based on the set of N — 1 angular variables defined in 
Eq.(2.11) 

a = («i(q), a 2 (q), • • • , a^-i(q)), (Al) 
each point q in the N- dimensional space will now be designated by 

(S, a 1 ,a 2 , ■ ■ -,q;jv-i), (A. 2) 

instead of (qi, q 2 , (fe, • • • , <7jv)- The corresponding line element can be written as 

d ~q=S h s dS + «j hjdaf, (A3) 
j=i 

the gradient is given by 

„ A 1 d ^ A 1 9 . . 



and 

T — - 

2 



T = - V (A5) 



can be decomposed into two parts: 

T = T S + T a , (A6) 

with 

1 d h a d 

5 = _ 2M^as ( 7^as ) ' ( ' 

m L vAf^Ai (ar\ 

2h s h a da/ h] da/' { ' 

N-l 

h a =n hj, (A9) 

i=i 

fc| = [(VS) 2 ]" 1 , fc? = [(Va 1 ) 2 ]- 1 ,"-, ^[(Va,-) 2 ]- 1 ,-. (A10) 
The volume element in the q-space is 

d N q = hsh a dSda (All) 

with 

JV-1 

da = J| day (A 12) 



in which 



and 
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Appendix B 



In the calculation of A n and r n in Section 3 the following expressions introduced in 
ref.[2] are needed. With the definition of C and D introduced in Eqs.(2.39) and (2.47), 
for n > 0, 

Cx 2n = [ S ^-x 2n = -^—x 2n . (5.1) 
Jo g2S g(2n) y ' 

n 

D[x 2n -r ln ] = (i + CT)- l c\x 2n - r, n ] = ^r m nX 2m . (5.2) 

m=l 

where 

1 r 2n-l 

5 J- n— 1 n 



#2n' 2g 2 (2n-2) 

\ — J 
to(2^T 



(2n-l)(2n-3)---(2m+l) r , „ x 

r «- = - ^7=T for l<"»<n-l (5.3) 



and in particular 

(2n-l)!! 

For n > 0, 

Or 2 " +1 = / 5 ^x 2n+1 = — -x 2n+l . (5.4) 

Jo g2S g(2n +1) v 1 

It is straightforward to establish the following: 

n 

Dx 2n+1 = (1 + CT)- l Cx 2n+1 = lm nx 2m+l (5.5) 

where 



m=0 



1 n 



In n / . -, \ j 7n— 1 n 



<?(2n + l)' ' n_in (? 2 (2n-l)' 
n — 1) • • • (to + 

# n - m+1 (2TO + 1) 



n(n- 1)---(to + 1) 
7m n = — rr^- for < to < n (5.6) 



and in particular 

n 

10n = 



g n+l 



It is convenient to extend the definitions of 7 m n and T m n , by defining 

7™ n = T m „ = for to > n (5.7) 

and 

r n = o. 
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with 



Appendix C 

For the special case of spherical coordinate in three-dimension and for the coulomb 
potential we have S = g 2 r. Correspondingly, 

a = (0, <p) (C.l) 

hs — 1, h& = r, h v — rsim? (C.2) 
For any function /(r) we have 

< W = {%? + + (0.3) 

T = T s + T a (C.4) 

_ 1 1 d d 
T s = r — , 

Ta = ^ (C.5) 

and L 2 is the operator of the square of the angular momentum. The Green functions used 
in the iteration are defined as 

G= (1 +DT a )~ 1 D = (1 + CT)~ 1 C. (C.6) 

Assuming \P(0) = $(0) and \I/(r — > oo) = 0, we have 

C = 1 / rfr (C.7) 

r oo 

D = -2[ dr^e^J rfrW 2 ^'. (C.8) 

r r 

In the iteration the following formula are applied: 

i n+1 

(1 + CT)- 1 ^ = - E V«r m + Ai,n(l + CT) _1 C • 1, (C.9) 

# m=2 
1 n+1 

Dr ' n = — A m,nr m + A ltn D ■ 1, (CIO) 

# m=2 
. _ 1 



where 



n+1 
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_ 1 (w + 2)(w + l)---(m + 2) 

A m ,n ~ (2 ^ 2)n+1 _ m m 

A ^ = j^(n + 2)(n + l)---3. (C7.ll) 

Although D ■ 1 = (1 + CT) _1 C ■ 1 is divergent, all terms including this factor are cancelled 
in the iteration. Defining cost? = £ we have 

r 2 T a C = £ 

and 

rXr = m(^±V - ^r- 2 )- (C.12) 
For Legendre function Pj(£), we have 

r 2 T a P,(0 = P(0- (C?.13) 

Based on the above formula the two iteration procedures could be performed. 
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